/*_glsl_gen.h - GLSL 'vec234' and 'mat234' datatypes in C++
 *
 *
 * Copyright (C) 2007 Bert van der Weerd <hypercharge@gmail.com>
 *
 * This software is licenced under the 'New BSD Licence' and part of
 * the iondrive project on google code.
 *
 * WARNING: Code status = Beta
 */


/*
 * --------
 * OVERVIEW
 * 
 * This code does not depend on or use OpenGL nor GLSL in any way. it
 * just uses the GLSL API.
 *
 * Because the vec<T,N> and mat<T,N> types are templated to both
 * element type T and dimension N, the following is also possible:
 *
 *   vec< complex<float>, 200 >  myVector;
 *
 * ---------
 * BEHAVIOUR
 *
 * In pseudo-BNF-code, the 'vec' and 'mat' types can do this:
 * 
 * 1) Create
 *	Matrices or Vectors
 * of type
 *	[bool|int|float|double|complex<T>]
 * of dimension
 *	[N].
 *
 * 2) Provide GLSL api as [global|member] functions, where 'GLOBAL API' is the
 * original
 * 
 *	A = transpose(B);
 *	
 * syntax, and 'MEMBER API' is the
 * 
 *	A = B;
 *	A.transpose();
 *
 * syntax. The latter syntax can be used to avoid temporary object creation,
 * which can be too costly for bigger size (vectors|matrices). You can
 * mix both syntax forms.
 *
 * 3) Provide [(Read|Write) (Binary|Text)] input/output of the two
 * types. (Only text-ostream is implemented)
 *	
 * 4) These classes are written with these main goals in mind:
 *
 *	1) Easy API that uses the well-known GLSL syntax.
 *	2) Speed and performance by using CLEAN CODE.
 *	3) performance is optimized for vec3/mat3 type of user code.
 *	3) for bigger N, optimization and scaling are difficult.
 *
 * ------------
 * CODE EXAMPLE
 * 
 *	vec3 x = init_vec(1.0f,2.0f,3.0f), y;
 *	mat3 A,B,C;
 *
 *	C = diag(x);
 *	A = transpose(C);
 *	B = A + C;
 *
 *	y = (A * B * C) * x;
 *
 *	return length(y);
 *
 * -----------------------------
 * ADDITIONAL NON-GLSL FUNCTIONS
 * 
 *   diag(vec<T,N> v) - create diagonal matrix from vector
 *   diag(T x, int N) - create diagonal (NxN)-matrix from scalair
 *   transpose(vec<T,N,M>) - I dont see (yet) why this one is missing
 *                           from the original glsl spec.
 *
 * ------
 * ISSUES
 * 
 * There are things in the GL Shading Language that cannot be
 * expressed in C++. These are the Swizzle Operations and Component
 * Name Aliasses. Just Remember: This is C++, not GLSL.
 *
 * The only functions NOT implemented are the texture functions, the
 * fragment shader specific functions, and the noise functions.
 *
 * In gcc, the not() builtin cannot be used. There is an operator!()
 * that performs the same function and is available everywhere. Also,
 * the memberwise constructors are moved into init_mat() and init_vec()
 * functions because constructors cannot be partially instantiated.
 *
 * MATRICES ARE IMPLEMENTED IN ROW-MAJOR ORDER, as a single
 * array, which is the standard in OpenGL.
 *
 * --------------------
 * IMPLEMENTATION NOTES
 *
 * The code #includes only <cmath> and only calls it explicitly from
 * the std namespace, like in 'std::sin(M_PI)'. The #include <iostream>
 * is just for simple mat/vec text output.
 *
 * The code consists of the following main classes:
 *
 *	template <typename T, int N> class vec { ... };
 *      template <typename T, int N, int M=N> class mat { ... };
 *
 * This allows us to write:
 *
 *      typedef vec<float,3> vec3;
 *      typedef mat<float,4> mat4;
 *
 * ------------------
 * OPTIMIZATION NOTES
 *
 * If you go for the last byte and cycle, you'll object to the glsl
 * syntax that it creates temporary copies, for example in the function
 * below:
 * 
 *      inline vec<T,N> cos(const vec<T,N>& v)
 *      { vec<T,N> out; for (int i=0;i<N;i++) out[i] = cos(v[i]); return out; }
 *
 * Now, for those who want to avoid this, the whole glsl api is also
 * accesable as memberfunctions into vec<T,N> and mat<T,N>, like this:
 * 
 *      inline void cos(const vec<T,N>& v)
 *      { for (int i=0;i<N;i++) (*this)[i] = cos(v[i]); }
 *
 * -------------
 * WEB RESOURCES
 *
 * * Matrix Layouts -- http://www.mindcontrol.org/~hplus/graphics/matrix-layout.html
 * * C++ Custom Stream I/O -- http://www.parashift.com/c++-faq-lite/input-output.html
 *
 * ----
 * TODO
 * 
 * 1) bring mat<> up to par with saru imp
 * 2) compare ColorRGB[A] with glsl api and bring up to par
 * 3) Read & up-to-par with http://www.morrowland.com/apron/tutorials/gl/gl_matrix.php
 *
 * -----
 * IDEAS
 * 
 * - stream reading of vec/mat
 * - combine stream read/write with ColMaj matrices so you can
 * - boolean vectors and logic for matrices and optimize using saru
 *   principles.
 * - implement Column Major matrices by partial templ. specialization.
 * - try out these classes a pure .NET code by vs2005 and check with
 *   saru for differences.
 * - std::vector integration might be a good idea? (low-priotiy).
 *
 */

#ifndef INCLUDED_GLSL_H
#define INCLUDED_GLSL_H

#ifndef __cplusplus
# error Who lives in a pineapple under the sea?
#endif


#include <iostream> // so we can read and write to streams
#include <cmath>

namespace glsl
{
	
// we still need to define M_PI and friends
#ifndef M_PI
# define M_PI 3.14159265358979323846
#endif
#ifndef M_LN2
# define M_LN2 0.693147180559945309417
#endif

// GNU C #defines log2, we dont want that
#ifdef __GNUC__
# ifdef log2
#  undef log2
# endif
#endif

// if min and max are #defines, remove them, we use our own
#ifdef min
# undef min
#endif
#ifdef max
# undef max
#endif




/************************************************************
 *
 * scalar builtin functions
 * 
 */


	//
	// Wrappers to standard math functions for type closure (these
	// get optimized away). They could also be used for faster
	// sin() etc functions by using tables.
	// 


	namespace cmath
	{
		template <typename T> inline T abs(const T& x)              { return std::abs(x); }
		template <typename T> inline T acos(const T& x)             { return std::acos(x); }
		template <typename T> inline T asin(const T& x)             { return std::asin(x); }
		template <typename T> inline T atan(const T& x)             { return std::atan(x); }
		template <typename T> inline T atan2(const T& x,const T& y) { return std::atan2(x,y); }
		template <typename T> inline T ceil(const T& x)             { return std::ceil(x); }
		template <typename T> inline T cos(const T& x)              { return std::cos(x); }
		template <typename T> inline T exp(const T& x)              { return std::exp(x); }
		template <typename T> inline T fabs(const T& x)             { return std::fabs(x); }
		template <typename T> inline T floor(const T& x)            { return std::floor(x); }
		template <typename T> inline T fmod(const T& x,const T& y)  { return std::fmod(x,y); }
		template <typename T> inline T log(const T& x)              { return std::log(x); }
		template <typename T> inline T pow(const T& x,const T& y)   { return std::pow(x,y); }
		template <typename T> inline T sin(const T& x)              { return std::sin(x); }
		template <typename T> inline T sqrt(const T& x)             { return std::sqrt(x); }
		template <typename T> inline T tan(const T& x)              { return std::tan(x); }
	}




	//
	// Scalair functions defined by GLSL
	// 


	template <typename T> inline T radians(const T& degrees)                            { return T(M_PI/180)*degrees; }
	template <typename T> inline T degrees(const T& radians)                            { return T(180/M_PI)*radians; }
	template <typename T> inline T exp2(const T& x)                                     { return T(cmath::exp(x * M_LN2)); }
	template <typename T> inline T log2(const T& x)                                     { return T(cmath::log(x) / M_LN2); }
	template <typename T> inline T inversesqrt(const T& x)                              { return 1/cmath::sqrt(x); }
	template <typename T> inline T sign(const T& x)                                     { return T((x>0) ? T(1) : ((x<0) ? T(-1):T(0))); }
	template <typename T> inline T fract(const T& x)                                    { return x - cmath::floor(x); }
	template <typename T> inline T mod(const T& x, const T& y)                          { return T(cmath::fmod(x, y)); }
	template <typename T> inline T min(const T& x, const T& y)                          { return y < x ? y : x; }
	template <typename T> inline T max(const T& x, const T& y)                          { return x < y ? y : x; }
	template <typename T> inline T clamp(const T& x, const T& minVal, const T& maxVal)  { return glsl::min(glsl::max(x,minVal),maxVal); }
	template <typename T> inline T mix(const T& x, const T& y, const T& a)              { return x*(1-a) + y*a; }
	template <typename T> inline T step(const T& edge, const T& x)                      { return x<=edge ? T(0) : T(1); }
	template <typename T> inline T smoothstep(const T& edge0,const T& edge1,const T& x) { T t = clamp((x-edge0) / (edge1-edge0), T(0), T(1)); return t * t * (3 - 2*t); }
	template <typename T> inline T length(const T& x)                                   { return cmath::sqrt(x*x); }
	template <typename T> inline T distance(const T& p0, const T& p1)                   { return length(p0-p1); }
	template <typename T> inline T dot(const T& x, const T& y)                          { return x*y; }
	template <typename T> inline T normalize(const T&)                                  { return T(1); } // this is not the most useful function in the world
	template <typename T> inline T faceforward(const T& N, const T& l, const T& Nref)   { return Nref*l < 0 ? N : -N; }
	template <typename T> inline T reflect(const T& I, const T& N)                      { return I - T(2)*N*I*N; }



/************************************************************
 *
 * vector class
 * 
 */



template <typename T, int N>
   class vec
{
	T data[N];
		
	public:

		inline vec()                         { for (int i=0;i<N;i++)   (*this)[i] = T(); }
		inline vec(const T& t)               { for (int i=0;i<N;i++)   (*this)[i] = t; }
		inline vec(vec<T,N-1> a, const T& t) { for (int i=0;i<N-1;i++) (*this)[i] = a[i];  (*this)[N-1] = t; } // dimension up cast
		inline vec(vec<T,N+1> a)             { for (int i=0;i<N;i++)   (*this)[i] = a[i];  }  // dimension down cast

		// *all* access is done through this operator
		inline T& operator[](const int i)       { return data[i]; }
		inline const T&  operator[](const int i) const { return data[i]; }

		
		
		//
		// experimental : non-temp-copy-creating interface
		//
		//
		// these inline functions below contain the glsl
		// functionality in member function form, which allows
		// you to use the vector class without generating the
		// temporary copies that inevitably arise when using
		// the glsl functional expression interface.
		//
		// ------------
		// CODE EXAMPLE
		//
		// vec3 x,y,z;
		// mat3 A,B;
		//
		// x = sin(A*y) + cos(B*z);
		//
		// ==[ BECOMES ]=======================
		// vec3 tmp1,tmp2,tmp3;
		//
		// tmp1.mul(A,y);                      // not implemented yet :(
		// tmp2.sin(tmp1);
		// tmp1.mul(B,z);
		// tmp3.cos(tmp1);
		// x.add(tmp2,tmp3);
		//
		//
		// 

		
		
		// init

		inline void init_vec(const T arg) { for (int i=0;i<N;i++) (*this)[i] = arg; }
		inline void init_vec(const T arg[]) { for (int i=0;i<N;i++) (*this)[i] = arg[i]; }
		inline void init_vec(const vec<T,N>& arg) { for (int i=0;i<N;i++) (*this)[i] = arg[i]; }

		
		// operators (result,v,a,b are vectors, t is scalair)

		
		// result.xxx(v)
		inline void add(const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] += rhs[i]; }
		inline void sub(const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] -= rhs[i]; }
		inline void mul(const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] *= rhs[i]; }
		inline void div(const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] /= rhs[i]; }
		// result.xxx(a,b);
		inline void add(const vec<T,N>& lhs,const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] + rhs[i]; }
		inline void sub(const vec<T,N>& lhs,const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] - rhs[i]; }
		inline void mul(const vec<T,N>& lhs,const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] * rhs[i]; }
		inline void div(const vec<T,N>& lhs,const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] / rhs[i]; }
		// result.xxx(t);
		inline void add(const T& rhs) { for (int i=0;i<N;i++) (*this)[i] += rhs; }
		inline void sub(const T& rhs) { for (int i=0;i<N;i++) (*this)[i] -= rhs; }
		inline void mul(const T& rhs) { for (int i=0;i<N;i++) (*this)[i] *= rhs; }
		inline void div(const T& rhs) { for (int i=0;i<N;i++) (*this)[i] /= rhs; }
		// result.xxx(v,t);
		inline void add(const vec<T,N>& lhs,const T& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] + rhs; }
		inline void sub(const vec<T,N>& lhs,const T& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] - rhs; }
		inline void mul(const vec<T,N>& lhs,const T& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] * rhs; }
		inline void div(const vec<T,N>& lhs,const T& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] / rhs; }

		inline void neg() { for (int i=0;i<N;i++) (*this)[i] = - ((*this)[i]); } // result.neg();
		inline void neg(const vec<T,N>& v) { for (int i=0;i<N;i++) (*this)[i] = - (v[i]); } // result.neg(vec_to_negate);


		// glsl functions
		

		inline void radians()                                                                   { for (int i=0;i<N;i++) (*this)[i] = radians((*this)[i]); }
		inline void radians(const vec<T,N>& v)                                                  { for (int i=0;i<N;i++) (*this)[i] = radians(v[i]); }
 		inline void degrees()                                                                   { for (int i=0;i<N;i++) (*this)[i] = degrees((*this)[i]); }
		inline void degrees(const vec<T,N>& v)                                                  { for (int i=0;i<N;i++) (*this)[i] = degrees(v[i]); }
		inline void sin()                                                                       { for (int i=0;i<N;i++) (*this)[i] = cmath::sin((*this)[i]); }
		inline void sin(const vec<T,N>& v)                                                      { for (int i=0;i<N;i++) (*this)[i] = cmath::sin(v[i]); }
		inline void cos()                                                                       { for (int i=0;i<N;i++) (*this)[i] = cmath::cos((*this)[i]); }
		inline void cos(const vec<T,N>& v)                                                      { for (int i=0;i<N;i++) (*this)[i] = cmath::cos(v[i]); }
		inline void tan()                                                                       { for (int i=0;i<N;i++) (*this)[i] = cmath::tan((*this)[i]); }
		inline void tan(const vec<T,N>& v)                                                      { for (int i=0;i<N;i++) (*this)[i] = cmath::tan(v[i]); }
		inline void asin()                                                                      { for (int i=0;i<N;i++) (*this)[i] = cmath::asin((*this)[i]); }
		inline void asin(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = cmath::asin(v[i]); }
		inline void acos()                                                                      { for (int i=0;i<N;i++) (*this)[i] = cmath::acos((*this)[i]); }
		inline void acos(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = cmath::acos(v[i]); }
		inline void atan()                                                                      { for (int i=0;i<N;i++) (*this)[i] = cmath::atan((*this)[i]); }
		inline void atan(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = cmath::atan(v[i]); }
		inline void atan2(const vec<T,N>& rhs)                                                  { for (int i=0;i<N;i++) (*this)[i] = cmath::atan2((*this)[i],rhs[i]); }
		inline void atan(const vec<T,N>& v,const vec<T,N>& u)                                   { for (int i=0;i<N;i++) (*this)[i] = cmath::atan2(v[i],u[i]); }
		inline void pow(const vec<T,N>& rhs)                                                    { for (int i=0;i<N;i++) (*this)[i] = cmath::pow((*this)[i],rhs[i]); }
		inline void pow(const vec<T,N>& lhs,const vec<T,N> rhs)                                 { for (int i=0;i<N;i++) (*this)[i] = cmath::pow(lhs[i],rhs[i]); }
		inline void exp2()                                                                      { for (int i=0;i<N;i++) (*this)[i] = exp2((*this)[i]); }
		inline void exp2(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = exp2(v[i]); }
		inline void log2()                                                                      { for (int i=0;i<N;i++) (*this)[i] = log2((*this)[i]); }
		inline void log2(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = log2(v[i]); }
		inline void exp()                                                                       { for (int i=0;i<N;i++) (*this)[i] = cmath::exp((*this)[i]); }
		inline void exp(const vec<T,N>& v)                                                      { for (int i=0;i<N;i++) (*this)[i] = cmath::exp(v[i]); }
		inline void log()                                                                       { for (int i=0;i<N;i++) (*this)[i] = cmath::log((*this)[i]); }
		inline void log(const vec<T,N>& v)                                                      { for (int i=0;i<N;i++) (*this)[i] = cmath::log(v[i]); }
		inline void sqrt()                                                                      { for (int i=0;i<N;i++) (*this)[i] = cmath::sqrt((*this)[i]); }
		inline void sqrt(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = cmath::sqrt(v[i]); }
		inline void inversesqrt()                                                               { for (int i=0;i<N;i++) (*this)[i] = inversesqrt((*this)[i]); }
		inline void inversesqrt(const vec<T,N>& v)                                              { for (int i=0;i<N;i++) (*this)[i] = inversesqrt(v[i]); }
		inline void abs()                                                                       { for (int i=0;i<N;i++) (*this)[i] = cmath::abs((*this)[i]); }
		inline void abs(const vec<T,N>& v)                                                      { for (int i=0;i<N;i++) (*this)[i] = cmath::abs(v[i]); }
		inline void sign()                                                                      { for (int i=0;i<N;i++) (*this)[i] = sign((*this)[i]); }
		inline void sign(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = sign(v[i]); }
		inline void floor()                                                                     { for (int i=0;i<N;i++) (*this)[i] = cmath::floor((*this)[i]); }
		inline void floor(const vec<T,N>& v)                                                    { for (int i=0;i<N;i++) (*this)[i] = cmath::floor(v[i]); }
		inline void ceil()                                                                      { for (int i=0;i<N;i++) (*this)[i] = cmath::ceil((*this)[i]); }
		inline void ceil(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = cmath::ceil(v[i]); }
		inline void fract()                                                                     { for (int i=0;i<N;i++) (*this)[i] = fract((*this)[i]); }
		inline void fract(const vec<T,N>& v)                                                    { for (int i=0;i<N;i++) (*this)[i] = fract(v[i]); }
		inline void mod(const T& rhs)                                                           { for (int i=0;i<N;i++) (*this)[i] = glsl::mod((*this)[i],rhs); }
		inline void mod(const vec<T,N>& lhs,const vec<T,N>& rhs)                                { for (int i=0;i<N;i++) (*this)[i] = glsl::mod(lhs[i],rhs[i]); }
		inline void min(const T& rhs)                                                           { for (int i=0;i<N;i++) (*this)[i] = glsl::min((*this)[i],rhs); }
		inline void min(const vec<T,N>& lhs,const vec<T,N>& rhs)                                { for (int i=0;i<N;i++) (*this)[i] = glsl::min(lhs[i],rhs[i]); }
		inline void max(const T& rhs)                                                           { for (int i=0;i<N;i++) (*this)[i] = glsl::max((*this)[i],rhs); }
		inline void max(const vec<T,N>& lhs,const vec<T,N>& rhs)                                { for (int i=0;i<N;i++) (*this)[i] = glsl::max(lhs[i],rhs[i]); }
		inline void clamp(const T& minVal, const T& maxVal)                                     { for (int i=0;i<N;i++) (*this)[i] = clamp((*this)[i],minVal,maxVal); }
		inline void clamp(const vec<T,N>& v, const T& minVal, const T& maxVal)                  { for (int i=0;i<N;i++) (*this)[i] = clamp(v[i],minVal,maxVal); }
		inline void clamp(const vec<T,N>& minVal, const vec<T,N>& maxVal)                       { for (int i=0;i<N;i++) (*this)[i] = clamp((*this)[i],minVal[i],maxVal[i]); }
		inline void clamp(const vec<T,N>& v,const vec<T,N>& minVal, const vec<T,N>& maxVal)     { for (int i=0;i<N;i++) (*this)[i] = clamp(v[i],minVal[i],maxVal[i]); }
		inline void mix(const vec<T,N>& a, const vec<T,N>& b, const T& t)                       { for (int i=0;i<N;i++) (*this)[i] = mix(a[i],b[i],t); }
		inline void mix(const vec<T,N>& a, const vec<T,N>& b, const vec<T,N>& c)                { for (int i=0;i<N;i++) (*this)[i] = mix(a[i],b[i],c[i]); }
		inline void step(const vec<T,N>& edge, const vec<T,N>& x)                               { for (int i=0;i<N;i++) (*this)[i] = step(edge[i],x[i]); }
		inline void smoothstep(const vec<T,N>& edge0, const vec<T,N>& edge1, const vec<T,N>& x) { for (int i=0;i<N;i++) (*this)[i] = smoothstep(edge0[i],edge1[i],x[i]); }

		inline void normalize()                  { T len =   length(); if (len == 0) return; this->div(len); }
		inline void normalize(const vec<T,N>& v) { T len = v.length(); if (len == 0) this->init(0); else this->div(v,len); }

		inline void faceforward(const vec<T,N>& n, const vec<T,N>& i, const vec<T,N>& nref)
		{
			// return dot(nref,i) < 0 ? n : -n;

			T tmp1 = nref.dot(i);
			if (tmp1 < 0)
				for(int tmp2=0;tmp2<N;tmp2++) (*this)[tmp2] = n[tmp2];
			else
				for(int tmp3=0;tmp3<N;tmp3++) (*this)[tmp3] = -n[tmp3];
		}
		inline void reflect(const vec<T,N>& i, const vec<T,N>& n)
		{
			// return i - 2*dot(n,i)*n;
			vec<T,N> tmp2;
			tmp2.mul(n,2*n.dot(i));
			this->sub(i,tmp2);
		}

		inline T length(const vec<T,N>& v)
		{
			return T(sqrt(dot(v,v)));
		}
		inline T distance(const vec<T,N>& a, const vec<T,N>& b)
		{
			// return length(a-b);
			vec<T,N> tmp1;
			tmp1.sub(a,b);
			return tmp1.length();
		}
		inline T dot(const vec<T,N>& rhs)
		{
			// T t; for (int i = 0; i < N; i++) t += a[i]*b[i]; return t;

			T t;
			for (int i=0;i<N;i++) t += (*this)[i]*rhs[i];
			return t;
		}
		inline void cross(const vec<T,3>& a,const vec<T,3>& b)
		{
			(*this)[0] = a[1]*b[2] - a[2]*b[1];
			(*this)[1] = a[2]*b[0] - a[0]*b[2];
			(*this)[2] = a[0]*b[1] - a[1]*b[0];
		}
		inline void cross(const vec<T,3>& b) // above cross() could be more efficient
		{
			vec<T,3> a;
			
			a.init_vec((*this)[0], (*this)[1], (*this)[2]);

			(*this)[0] = a[1]*b[2] - a[2]*b[1];
			(*this)[1] = a[2]*b[0] - a[0]*b[2];
			(*this)[2] = a[0]*b[1] - a[1]*b[0];
		}

		// Additional saru functions
		inline T squared_length() { T tmp; for (int i=0;i<N;i++) tmp += (*this)[i]*(*this)[i]; return tmp; }
		inline void unit_cross(const vec<T,3>& b)
		{
			vec<T,3> result;

			result[0] = (*this)[1]*b[2] - (*this)[2]*b[1];
			result[1] = (*this)[2]*b[0] - (*this)[0]*b[2];
			result[2] = (*this)[0]*b[1] - (*this)[1]*b[0];

			result.normalize();

			return result;
		}
};


// C++ LIMIT: We cannot do partial template specialization in the
// constructor, so these init_vec[234]() functions are provided as a fix.

template <typename T> inline vec<T,2> init_vec(const T& x, const T& y)                         { vec<T,2> out; out[0] = x; out[1] = y; return out; }
template <typename T> inline vec<T,3> init_vec(const T& x, const T& y, const T& z)             { vec<T,3> out; out[0] = x; out[1] = y; out[2] = z; return out; }
template <typename T> inline vec<T,4> init_vec(const T& x, const T& y, const T& z, const T& w) { vec<T,4> out; out[0] = x; out[1] = y; out[2] = z; out[3] = w; return out; }

// REQUIRED: Vector operators

template <typename T, int N> inline vec<T,N> operator+(const vec<T,N>& a, const vec<T,N>& b) { vec<T,N> out; for (int i=0;i<N;i++) out[i] = a[i]+b[i]; return out; }
template <typename T, int N> inline vec<T,N> operator-(const vec<T,N>& a, const vec<T,N>& b) { vec<T,N> out; for (int i=0;i<N;i++) out[i] = a[i]-b[i]; return out; }
template <typename T, int N> inline vec<T,N> operator*(const vec<T,N>& a, const vec<T,N>& b) { vec<T,N> out; for (int i=0;i<N;i++) out[i] = a[i]*b[i]; return out; }
template <typename T, int N> inline vec<T,N> operator/(const vec<T,N>& a, const vec<T,N>& b) { vec<T,N> out; for (int i=0;i<N;i++) out[i] = a[i]/b[i]; return out; }

template <typename T, int N> inline vec<T,N> operator*(const T& t, const vec<T,N>& v)        { vec<T,N> out; for (int i=0;i<N;i++) out[i] = t*v[i]; return out; }
template <typename T, int N> inline vec<T,N> operator*(const vec<T,N>& v, const T& t)        { vec<T,N> out; for (int i=0;i<N;i++) out[i] = v[i]*t; return out; }
template <typename T, int N> inline vec<T,N> operator+(const T& t, const vec<T,N>& v)        { vec<T,N> out; for (int i=0;i<N;i++) out[i] = t+v[i]; return out; }
template <typename T, int N> inline vec<T,N> operator+(const vec<T,N>& v, const T& t)        { vec<T,N> out; for (int i=0;i<N;i++) out[i] = v[i]+t; return out; }
template <typename T, int N> inline vec<T,N> operator-(const T& t, const vec<T,N>& v)        { vec<T,N> out; for (int i=0;i<N;i++) out[i] = t-v[i]; return out; }
template <typename T, int N> inline vec<T,N> operator-(const vec<T,N>& v, const T& t)        { vec<T,N> out; for (int i=0;i<N;i++) out[i] = v[i]-t; return out; }
template <typename T, int N> inline vec<T,N> operator/(const vec<T,N>& v, const T& t)        { vec<T,N> out; T tmp=1/t; for (int i=0;i<N;i++) out[i] = v[i]*tmp; return out; } // speedup from saru

template <typename T, int N> inline vec<T,N>& operator+=(vec<T,N>& a, const vec<T,N>& b)     { a = a + b; return a; } // this looks improvable
template <typename T, int N> inline vec<T,N>& operator-=(vec<T,N>& a, const vec<T,N>& b)     { a = a - b; return a; }
template <typename T, int N> inline vec<T,N>& operator*=(vec<T,N>& a, const vec<T,N>& b)     { a = a * b; return a; }
template <typename T, int N> inline vec<T,N>& operator/=(vec<T,N>& a, const vec<T,N>& b)     { a = a / b; return a; }
template <typename T, int N> inline vec<T,N>& operator+=(vec<T,N>& a, const T& t)            { a = a + t; return a; }
template <typename T, int N> inline vec<T,N>& operator-=(vec<T,N>& a, const T& t)            { a = a - t; return a; }
template <typename T, int N> inline vec<T,N>& operator*=(vec<T,N>& a, const T& t)            { a = a * t; return a; }
template <typename T, int N> inline vec<T,N>& operator/=(vec<T,N>& a, const T& t)            { a = a / t; return a; }
// unary minus operator
template <typename T, int N> inline vec<T,N> operator-(const vec<T,N>& v)                    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = -v[i]; return out; }

// REQUIRED: Builtin GLSL functions

template <typename T, int N> inline vec<T,N> radians(const vec<T,N>& v)                                                  { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = radians(v[i]); return out; }
template <typename T, int N> inline vec<T,N> degrees(const vec<T,N>& v)                                                  { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = degrees(v[i]); return out; }
template <typename T, int N> inline vec<T,N> sin(const vec<T,N>& v)                                                      { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = cmath::sin(v[i]); return out; }
template <typename T, int N> inline vec<T,N> cos(const vec<T,N>& v)                                                      { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = cmath::cos(v[i]); return out; }
template <typename T, int N> inline vec<T,N> tan(const vec<T,N>& v)                                                      { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = cmath::tan(v[i]); return out; }
template <typename T, int N> inline vec<T,N> asin(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = cmath::asin(v[i]); return out; }
template <typename T, int N> inline vec<T,N> acos(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = cmath::acos(v[i]); return out; }
template <typename T, int N> inline vec<T,N> atan(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = cmath::atan(v[i]); return out; }
template <typename T, int N> inline vec<T,N> atan(const vec<T,N>& y,const vec<T,N>& x)                                   { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = cmath::atan2(y[i],x[i]); return out; }
template <typename T, int N> inline vec<T,N> pow(const vec<T,N>& a,const vec<T,N>& x)                                    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = cmath::pow(a[i],x[i]); return out; }
template <typename T, int N> inline vec<T,N> exp2(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = exp2(v[i]); return out; }
template <typename T, int N> inline vec<T,N> log2(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = log2(v[i]); return out; }
template <typename T, int N> inline vec<T,N> exp(const vec<T,N>& v)                                                      { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = cmath::exp(v[i]); return out; }
template <typename T, int N> inline vec<T,N> log(const vec<T,N>& v)                                                      { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = cmath::log(v[i]); return out; }
template <typename T, int N> inline vec<T,N> sqrt(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = cmath::sqrt(v[i]); return out; }
template <typename T, int N> inline vec<T,N> inversesqrt(const vec<T,N>& v)                                              { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = inversesqrt(v[i]); return out; }
template <typename T, int N> inline vec<T,N> abs(const vec<T,N>& v)                                                      { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = cmath::fabs(v[i]); return out; }
template <typename T, int N> inline vec<T,N> sign(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = sign(v[i]); return out; }
template <typename T, int N> inline vec<T,N> floor(const vec<T,N>& v)                                                    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = cmath::floor(v[i]); return out; }
template <typename T, int N> inline vec<T,N> ceil(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = cmath::ceil(v[i]); return out; }
template <typename T, int N> inline vec<T,N> fract(const vec<T,N>& v)                                                    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = fract(v[i]); return out; }
template <typename T, int N> inline vec<T,N> mod(const vec<T,N>& v, const T& t)                                          { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = glsl::mod(v[i],t); return out; }
template <typename T, int N> inline vec<T,N> mod(const vec<T,N>& y,const vec<T,N>& x)                                    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = glsl::mod(y[i],x[i]); return out; }
template <typename T, int N> inline vec<T,N> min(const vec<T,N>& v, const T& t)                                          { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = glsl::min(v[i],t); return out; }
template <typename T, int N> inline vec<T,N> min(const vec<T,N>& y,const vec<T,N>& x)                                    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = glsl::min(y[i],x[i]); return out; }
template <typename T, int N> inline vec<T,N> max(const vec<T,N>& v, const T& t)                                          { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = glsl::max(v[i],t); return out; }
template <typename T, int N> inline vec<T,N> max(const vec<T,N>& y,const vec<T,N>& x)                                    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = glsl::max(y[i],x[i]); return out; }
template <typename T, int N> inline vec<T,N> clamp(const vec<T,N>& v, const T& minVal, const T& maxVal)                  { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = clamp(v[i],minVal,maxVal); return out; }
template <typename T, int N> inline vec<T,N> clamp(const vec<T,N>& v, const vec<T,N>& minVal, const vec<T,N>& maxVal)    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = clamp(v[i],minVal[i],maxVal[i]); return out; }
template <typename T, int N> inline vec<T,N> mix(const vec<T,N>& a,const vec<T,N>& b, const T& t)                        { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = mix(a[i],b[i],t); return out; }
template <typename T, int N> inline vec<T,N> mix(const vec<T,N>& a,const vec<T,N>& b, const vec<T,N>& c)                 { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = mix(a[i],b[i],c[i]); return out; }
template <typename T, int N> inline vec<T,N> step(const vec<T,N>& edge, const vec<T,N>& x)                               { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = step(edge[i],x[i]); return out; }
template <typename T, int N> inline vec<T,N> smoothstep(const vec<T,N>& edge0, const vec<T,N>& edge1, const vec<T,N>& x) { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = smoothstep(edge0[i],edge1[i],x[i]); return out; }

template <typename T, int N> inline vec<T,N> normalize(const vec<T,N>& v)                                                { T len = length(v); if (len == 0) return v; return v/len; }
template <typename T, int N> inline vec<T,N> faceforward(const vec<T,N>& n, const vec<T,N>& i, const vec<T,N>& nref)     { return dot(nref,i) < 0 ? n : -n; }
template <typename T, int N> inline vec<T,N> reflect(const vec<T,N>& i, const vec<T,N>& n)                               { return i - dot(n,i)*n*T(2); }

template <typename T, int N> inline T length(const vec<T,N>& v)                                                          { return T((cmath::sqrt(dot(v,v)))); }
template <typename T, int N> inline T distance(const vec<T,N>& a, const vec<T,N>& b)                                     { return length(a-b); }
template <typename T, int N> inline T dot(const vec<T,N>& a, const vec<T,N>& b)                                          { T tmp; for (int i=0;i<N;i++) tmp += a[i]*b[i]; return tmp; }

// partial specialization for CrossProduct: it only exists in 3D vectors
template <typename T> inline vec<T,3> cross(const vec<T,3>& a,const vec<T,3>& b) {
	vec<T,3> result;
	
	result[0] = a[1]*b[2] - a[2]*b[1];
	result[1] = a[2]*b[0] - a[0]*b[2];
	result[2] = a[0]*b[1] - a[1]*b[0];

	return result;
}

/* boolean vectors - there are currently no plans to try to implement
 * these boolean vectors without the use of temporaries. If speed on
 * these functions becomes an issue we'll see about that :-)
 */

// partial specialization for boolean vectors
template <typename T, int N> inline vec<bool,N> lessThan(const vec<T,N>& a, const vec<T,N>& b)         { vec<bool,N> out; for (int i = 0; i < N; i++) out[i] = a[i] < b[i]; return out; }
template <typename T, int N> inline vec<bool,N> lessThanEqual(const vec<T,N>& a, const vec<T,N>& b)    { vec<bool,N> out; for (int i = 0; i < N; i++) out[i] = a[i] <= b[i]; return out; }
template <typename T, int N> inline vec<bool,N> greaterThan(const vec<T,N>& a, const vec<T,N>& b)      { vec<bool,N> out; for (int i = 0; i < N; i++) out[i] = a[i] > b[i]; return out; }
template <typename T, int N> inline vec<bool,N> greaterThanEqual(const vec<T,N>& a, const vec<T,N>& b) { vec<bool,N> out; for (int i = 0; i < N; i++) out[i] = a[i] >= b[i]; return out; }
template <typename T, int N> inline vec<bool,N> equal(const vec<T,N>& a, const vec<T,N>& b)            { vec<bool,N> out; for (int i = 0; i < N; i++) out[i] = a[i] == b[i]; return out; }
template <typename T, int N> inline vec<bool,N> notEqual(const vec<T,N>& a, const vec<T,N>& b)         { vec<bool,N> out; for (int i = 0; i < N; i++) out[i] = a[i] != b[i]; return out; }

template <int N> inline bool any(const vec<bool,N>& v) { for (int i=0;i<N;i++) if ( v[i]) return true;  return false; }
template <int N> inline bool all(const vec<bool,N>& v) { for (int i=0;i<N;i++) if (!v[i]) return false; return true;  }

// C++ LIMIT: Because gcc reserves not() we cannot use that GLSL func. we
// use operator!() instead.
#ifndef __GNUC__
template <int N> inline vec<bool,N> not(const vec<bool,N>& v) {
	vec<bool,N> out;
	for (int i = 0; i < N; i++) out[i] = !v[i];
	return out;
}
#endif
// this operator!() can be used instead of the above not()
template <int N> inline vec<bool,N> operator!(const vec<bool,N>& v) {
	vec<bool,N> out;
	for (int i = 0; i < N; i++) out[i] = !v[i];
	return out;
}

//
// additional non-standard functions
//

template <typename T,int N> inline T squared_length(const vec<T,N>& v)   { T tmp; for (int i=0;i<N;i++) tmp += v[i]*v[i]; return tmp; }
template <typename T> inline vec<T,3> unit_cross(const vec<T,3>& a,const vec<T,3>& b)
{
	vec<T,3> result;

	result[0] = a[1]*b[2] - a[2]*b[1];
	result[1] = a[2]*b[0] - a[0]*b[2];
	result[2] = a[0]*b[1] - a[1]*b[0];

	result.normalize();
	
	return result;
}




/************************************************************
 *
 * matrix class
 * 
 */




template <typename T, int N, int M = N>
   class mat
{
	T data[N*M];

	public:

		inline mat() { for (int i = 0; i < N*M; i++) (*this)[i] = T(); }

		// *all* access is done through this operator
		inline T& operator[](const int i)       { return data[i]; }
		inline const T&  operator[](const int i) const { return data[i]; }

		inline T& element(const int i, const int j) { return (*this)[i*N+j]; }
		inline const T& element(const int i, const int j) const { return (*this)[i*N+j]; }
};
   
// use init() to initialize explicitly the 2,3,4D case

template <typename T> inline mat<T,2> init_mat(const T& m0,
					       const T& m1,
					       const T& m2,
					       const T& m3)
{
	mat<T,2> m;

	m[0] = m0;
	m[1] = m1;
	m[2] = m2;
	m[3] = m3;

	return m;
}
template <typename T> inline mat<T,3> init_mat(const T& m0,
					       const T& m1,
					       const T& m2,
					       const T& m3,
					       const T& m4,
					       const T& m5,
					       const T& m6,
					       const T& m7,
					       const T& m8)
{
	mat<T,3> m;

	m[0]=m0;
	m[1]=m1;
	m[2]=m2;
	m[3]=m3;
	m[4]=m4;
	m[5]=m5;
	m[6]=m6;
	m[7]=m7;
	m[8]=m8;

	return m;
}
template <typename T> inline mat<T,4> init_mat(const T& m0,
					       const T& m1,
					       const T& m2,
					       const T& m3,
					       const T& m4,
					       const T& m5,
					       const T& m6,
					       const T& m7,
					       const T& m8,
					       const T& m9,
					       const T& m10,
					       const T& m11,
					       const T& m12,
					       const T& m13,
					       const T& m14,
					       const T& m15)
{
	mat<T,4> m;

	m[0]=m0;
	m[1]=m1;
	m[2]=m2;
	m[3]=m3;
	m[4]=m4;
	m[5]=m5;
	m[6]=m6;
	m[7]=m7;
	m[8]=m8;
	m[9]=m9;
	m[10]=m10;
	m[11]=m11;
	m[12]=m12;
	m[13]=m13;
	m[14]=m14;
	m[15]=m15;

	return m;
}

// REQUIRED: matrix operations

template <typename T, int N> inline mat<T,N> operator+(mat<T,N> a, mat<T,N> b) { mat<T,N> out; for (int i = 0; i < N*N; i++) out[i] = a[i]+b[i]; return out; }
template <typename T, int N> inline mat<T,N> operator-(mat<T,N> a, mat<T,N> b) { mat<T,N> out; for (int i = 0; i < N*N; i++) out[i] = a[i]-b[i]; return out; }
template <typename T, int N> inline mat<T,N> operator+(const T& a, mat<T,N> m) { mat<T,N> out; for (int i = 0; i < N*N; i++) out[i] = a+m[i]; return out; }
template <typename T, int N> inline mat<T,N> operator+(mat<T,N> m, const T& a) { return a+m; }
template <typename T, int N> inline mat<T,N> operator-(const T& a, mat<T,N> m) { mat<T,N> out; for (int i = 0; i < N*N; i++) out[i] = a-m[i]; return out; }
template <typename T, int N> inline mat<T,N> operator-(mat<T,N> m, const T& a) { mat<T,N> out; for (int i = 0; i < N*N; i++) out[i] = m[i]-a; return out; }
template <typename T, int N> inline mat<T,N> operator*(const T& a, mat<T,N> m) { mat<T,N> out; for (int i = 0; i < N*N; i++) out[i] = a*m[i]; return out; }
template <typename T, int N> inline mat<T,N> operator*(mat<T,N> m, const T& a) { return a*m; }

// matrix-vector multiplication
template <typename T, int N> inline vec<T,N> operator*(mat<T,N> m, const vec<T,N>& v)
{
	vec<T,N> out;

	for (int j = 0; j < N; j++)
	{
		out[j] = T();
		
		for (int i = 0; i < N; i++)
			out[j] += m[j+N*i] * v[i];
	}
	
	return out;
}

// matrix-matrix multiplication
template <typename T, int N> inline mat<T,N> operator*(mat<T,N> a, mat<T,N> b)
{
	mat<T,N> out;

	for (int k = 0; k < N; k++)
	{
		for (int j = 0; j < N; j++)
		{
			out[k*N+j] = 0;
			
			for (int i = 0; i < N; i++)
				out[k*N+j] += a[j*N+i] * b[k*N+i];
		}
	}
	
	return out;
}

// builtins
template <typename T, int N, int M> inline mat<T,N,M> matrixcompmult(mat<T,N,M> a, mat<T,N,M> b) { mat<T,N,M> out; for (int i = 0; i < N*M; i++) out[i] = a[i]*b[i]; return out; }

//
// additional non-standard functions
//

template <typename T, int N, int M> inline mat<T,M,N> transpose(mat<T,N,M> m)
{
	mat<T,M,N> out;

	for (int j = 0; j < M; j++)
		for (int i = 0; i < N; i++)
			out[j*M+i] = m[i*N+j];
	
	return out;
}

template <typename T, int N> inline mat<T,N> diag(const vec<T,N>& v)
{
	mat<T,N> out;			// default ctor sets all to 0
	
	for (int i = 0; i < N; i++)
		out[i*N+i] = v[i];
	
	return out;
}
template <typename T, int N> inline mat<T,N> diag(const T& v, const int n=N)
{
	mat<T,N> out;			// mat3 foo = diag(1,3);=makeidenity; foo=diag(0,3); == makezero

	for (int i = 0; i < N; i++)
		out[i*N+i] = v[i];

	return out;
}




/************************************************************
 *
 * These are all the standard GLSL types as typedefs
 * 
 */




typedef vec<float,2>  vec2;
typedef vec<int,  2> ivec2;
typedef vec<bool, 2> bvec2;

typedef vec<float,3>  vec3;
typedef vec<int,  3> ivec3;
typedef vec<bool, 3> bvec3;

typedef vec<float,4>  vec4;
typedef vec<int,  4> ivec4;
typedef vec<bool, 4> bvec4;

typedef mat<float,2>  mat2;
typedef mat<float,3>  mat3;
typedef mat<float,4>  mat4;


/************************************************************
 *
 * The following part just handles text output for streams
 *
 */

template<typename T, int N> std::ostream& operator<<(std::ostream& os, vec<T,N> v)
{
	os << '(' << v[0];
	for (int i=1;i<N;i++) os << ',' << v[i];
	os << ')';

	return os;
}
template<typename T, int N, int M> std::ostream& operator<<(std::ostream& os, mat<T,N,M> A)
{
	os << '(';
	
	os << '(' << A.element(0,0);
	for (int i=1;i<M;i++) os << ',' << A.element(0,i);
	os << ')';
	
	for (int rows = 1; rows < N; rows++)
	{
		os << ',' << ' ';
		
		os << '(' << A.element(rows,0);
		for (int i=1;i<M;i++) os << ',' << A.element(rows,i);
		os << ')';
	}
	

	os << ')';
	
	return os;
}



} // end namespace glsl

#endif /* INCLUDED_GLSL_H */
